
/* 	This do-file replicates the results in the supplementary 
	appendix. 
	
	Please see README.txt for details on the data files. 
	
	Set the line "cd" below to the directory where the 
	data files have been extracted. 
	
	This do file requires an ado file -rdrobust- to run. 
	Type findit rdrobust to install. 
	
	IMPORTANT: The results were obtained with rdrobust version 6.0, 
		14Oct2014; a different version may give somewhat different 
		results. 
	
	All computations were done in Stata 12. 
	
	Please direct any questions, comments, or spotted errors to: 
	marko.klasnja@gmail.com
	
	*/

clear *
set more off
cd "C:/Users/`c(username)'/Dropbox/Romania/do files/disadv_JOP_replication/final"
set matsize 8000
set seed 7

* store results in a matrix
program define stores
	args y z /* y = matrix column z = matrix name */
	matrix `z'[1,`y'] = e(tau_cl)
	matrix `z'[2,`y'] = e(se_cl)
	matrix `z'[3,`y'] = e(pv_cl)
	matrix `z'[4,`y'] = e(h_bw)
	matrix `z'[5,`y'] = e(N)
	matrix `z'[6,`y'] = e(N_l)
	matrix `z'[7,`y'] = e(N_r)
	matrix rownames `z' = "Estimate" "SE" "p-value" "Bandwidth" "N" "N-below" "N-above"
end

* invert rows and columns when storing results
program define storeinv
	args y z /* y = matrix row, z = matrix name */
	matrix `z'[`y',1] = e(tau_cl)
	matrix `z'[`y',2] = e(se_cl)
	matrix `z'[`y',3] = e(pv_cl)
	matrix `z'[`y',4] = e(h_bw)
	matrix `z'[`y',5] = e(N)
	matrix `z'[`y',6] = e(N_l)
	matrix `z'[`y',7] = e(N_r)
	matrix colnames `z' = "Estimate" "SE" "p-value" "Bandwidth" "N" "N-below" "N-above"
end

* store resutls form McCrary tests
program define storemcc
	args x y z /* x = variable, y = matrix row, z = matrix name */
	local est = r(theta)
	local ster = r(se)
	local absest = abs(`est')
	local pval = (1-(normal(`absest'/`ster')))*2
	local bw = r(bandwidth)
	quietly sum `x' if `x' > -`bw' & `x' <= 0
	local N_l = r(N)
	quietly sum `x' if `x' > 0 & `x' <= `bw'
	local N_r = r(N)
	local Nall = `N_l' + `N_r'
	matrix `z'[`y',1] = `est'
	matrix `z'[`y',2] = `ster'
	matrix `z'[`y',3] = `pval'
	matrix `z'[`y',4] = `bw'
	matrix `z'[`y',5] = `Nall'
	matrix `z'[`y',6] = `N_l'
	matrix `z'[`y',7] = `N_r'
	matrix colnames `z' = "Estimate" "SE" "p-value" "Bandwidth" "N" "N-below" "N-above"
end

* store results for bandwidth sensitivity test across salary threshold
capture program drop stores_mw
program define stores_mw
	args y z /* y = matrix row, z = matrix name */
	matrix `z'[`y',1] = e(tau_cl)
	matrix `z'[`y',2] = e(se_cl)
	matrix `z'[`y',3] = e(N)
end

* store results from 2nd stage of fuzzy RD results
program define stores_f
	args y z /* y = matrix column z = matrix name */
	matrix `z'[1,`y'] = e(tau_F_cl)
	matrix `z'[2,`y'] = e(se_F_cl)
	matrix `z'[3,`y'] = e(pv_F_cl)
	matrix `z'[7,`y'] = e(h_bw)
	matrix `z'[8,`y'] = e(N)
	matrix `z'[9,`y'] = e(N_l)
	matrix `z'[10,`y'] = e(N_r)
end

* store results from Gelman-King models
program define storegk
	args x y /* x = row number, y = column number */
	matrix ptable = r(table)
	local j = `x' + 1
	local k = `x' + 2
	matrix gkmodels[`x',`y'] = _b[I2]
	matrix gkmodels[`j',`y'] = _se[I2]
	matrix gkmodels[`k',`y'] = ptable[4,3]
end


/*---------------------------------------------------------------------------*/
/* TABLE A1: ALTERNATIVE RUNNING VARIABLE */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

matrix res_alt = J(7,2,0)

* column 1: victory at t+1
qui rdrobust win alt_vote_margin
stores 1 res_alt

* column 2: vote margin at t+1
qui rdrobust win alt_vote_margin if largest == 1
stores 2 res_alt

* show results
matrix list res_alt

	
/*---------------------------------------------------------------------------*/
/* TABLE A2: FUZZY RDD MAIN ESTIMATES */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

matrix res_fuzzy = J(10,2,.)

* column 1: victory at t+1
rdrobust win vote_margin, fuzzy(win_t)
stores_f 1 res_fuzzy

* column 2: vote margin at t+1
rdrobust vote vote_margin, fuzzy(win_t)
stores_f 2 res_fuzzy

* show results
matrix list res_fuzzy


/*---------------------------------------------------------------------------*/
/* TABLE A3: GELMAN-KING (1990) CANDIDATE INCUMBENCY DISADVANTAGE */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

* build Gelman-King estimator
gen I2 = .
replace I2 = 1 if not_run_indiv == 0 & win_t == 1
replace I2 = -1 if not_run_indiv == 0 & win_t == 0
replace I2 = 0 if not_run_indiv == 1

gen P2 = 1 if win_t == 1
replace P2 = -1 if win_t == 0

* overall incumbency disadvantage in entire sample
reg vote svote1 P2 I2, r

* store results in a table
matrix gkmodels = J(12,2,.)

* panel 1, column 1: overall incumbency disadvantage within optimal bandwidth
* (using RDD to get the optimal bandwidth)
qui rdrobust vote vote_margin
global bw = e(h_bw)
qui reg vote svote1 P2 I2 if vote_margin > -$bw & vote_margin < $bw, r
storegk 2 1

* panel 1, column 2: overall incumbency disadvantage within smaller bandwidth
qui reg vote svote1 P2 I2 if vote_margin > -.1 & vote_margin < .1, r
storegk 2 2

/* Across salary threshold */

* panel 2, column 1: below salary threshold, within balanced pop. window
global pop_bw = .186
qui reg vote svote1 P2 I2 if pop_margin > -$pop_bw & pop_margin < 0 ///
	& vote_margin > -.1 & vote_margin < .1, r
storegk 6 1

* panel 2, column 2: above salary threshold, within balanced pop. window
qui reg vote svote1 P2 I2 if pop_margin > 0 & pop_margin < $pop_bw ///
	& vote_margin > -.1 & vote_margin < .1, r
storegk 6 2

/* Multiple-term vs. first-term incumbents */

* panel 3, column 1: first-term incumbents
qui reg vote svote1 P2 I2 if inc_prev == 0 & vote_margin > -.1 & vote_margin < .1
storegk 10 1

* panel 3, column 2: multiple-term incumbents
qui reg vote svote1 P2 I2 if inc_prev == 1 & vote_margin > -.1 & vote_margin < .1
storegk 10 2

matrix rownames gkmodels = "" "Estimate" "SE" "p-value" "" "Estimate" "SE" "p-value" ///
	"" "Estimate" "SE" "p-value"
matrix list gkmodels


/*---------------------------------------------------------------------------*/
/* FIGURE A2: 2008 MAYOR'S SALARY CHANGE ACROSS 7,000 THRESHOLD */
/*---------------------------------------------------------------------------*/

use 2008_salary_data, clear

* remove two outliers with very high salaries above the 
* threshold for easier inspection (result unaffected substantively)
rdbinselect salary pop_margin if city_code_unq ~= 1815 & city_code_unq ~= 2301, ///
	scale(6) p(3) lowerend(-.4) upperend(.4) ///
	graph_options(scheme(s1mono) title("") xtitle("Population") ///
	ytitle("Salary (Romanian Lei/Year)") legend(off) aspect(1) ///
	xlabel(-.4(.2).4) ylabel(, angle(0)) xline(0, lcolor(red))))
gr_edit xaxis1.major.num_rule_ticks = 5
gr_edit xaxis1.edit_tick 1 -0.4 `"4,200"', tickset(major)
gr_edit xaxis1.major.num_rule_ticks = 4
gr_edit xaxis1.edit_tick 1 -0.2 `"5,600"', tickset(major)
gr_edit xaxis1.major.num_rule_ticks = 3
gr_edit xaxis1.edit_tick 1 0 `"7000"', tickset(major)
gr_edit xaxis1.major.num_rule_ticks = 2
gr_edit xaxis1.edit_tick 5 0 `"7,000"', tickset(major)
gr_edit xaxis1.major.num_rule_ticks = 2
gr_edit xaxis1.edit_tick 1 0.2 `"8,400"', tickset(major)
gr_edit xaxis1.major.num_rule_ticks = 1
gr_edit xaxis1.edit_tick 1 0.4 `"9,800"', tickset(major)
graph export "$res_path/salary_jump.pdf", replace

* rd results
* raw data
rdrobust salary pop_margin

* logged salary
gen lnsalary = ln(salary)
rdrobust lnsalary pop_margin


/*---------------------------------------------------------------------------*/
/* FIGURE A3: SENSITIVITY TO BANDWIDTH SIZE. */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

matrix res_bw_size = J(28,4,0)

qui rdrobust win vote_margin
local bwmain = e(h_bw)

local count = 0
forval i = 3/30 {
	local k = `i' - 2
	local j = `i'/100
	qui rdrobust win vote_margin, h(`j')
	matrix res_bw_size[`k',1] = e(tau_cl)
	matrix res_bw_size[`k',2] = e(tau_cl) - 1.96*e(se_cl)
	matrix res_bw_size[`k',3] = e(tau_cl) + 1.96*e(se_cl)
	matrix res_bw_size[`k',4] = `j'
	}

svmat res_bw_size
twoway (line res_bw_size1 res_bw_size4, lcolor(black)) ///
	(line res_bw_size2 res_bw_size4, lpattern(dash) lcolor(black)) ///
	(line res_bw_size3 res_bw_size4, lpattern(dash) lcolor(black)), ///
	scheme(s1mono) ylabel(, nogrid) yline(0, lcolor(gs13)) legend(off) ///
	aspect(1) xtitle(Bandwidth size) ytitle(Incumbency Disadvantage Estimate) ///
	yscale(range(-.5(.1).5)) ylabel(#7) xline(`bwmain', lcolor(red))

	
/*---------------------------------------------------------------------------*/
/* TABLE A5: PRETREATMENT BALANCE TESTS -- VOTE MARGIN */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

matrix pretreat = J(13,7,.)

* row 1: vote margin at t-1
qui rdrobust vote_margin_t1 vote_margin
storeinv 1 pretreat

* row 2: vote margin at t-1 around population/salary cutoff
global salary_bw = .186
qui rdrobust vote_margin_t1 vote_margin if pop_margin > -$salary_bw & pop_margin < $salary_bw
storeinv 2 pretreat

* row 3: victory at t-2
qui rdrobust win_t2 vote_margin
storeinv 3 pretreat

* row 4: victory at t-2 around population/salary cutoff
qui rdrobust win_t2 vote_margin if pop_margin > -$salary_bw & pop_margin < $salary_bw
storeinv 4 pretreat

* rows 5-13: other electoral and municipal covariates (at t)
replace area_agr_t = area_agr_t/area_total_t
local vars turnout1_t effround1_t local_cop_t pop_t empl_t educ_t area_agr_t r_total e_total
local n_vars : word count `vars'

forval i = 1/`n_vars' {
	qui rdrobust `: word `i' of `vars'' vote_margin
	local j = `i' + 4
	storeinv `j' pretreat
	}
	
matrix rownames pretreat = vote_margin_t1 vote_margin_t1_pop win_t2 win_t2_pop ///
	turnout eff_num_par copartisan population empl schooling agr_area ///
	revenues expenditures

* show results
matrix list pretreat


/*---------------------------------------------------------------------------*/
/* TABLE A6: PRETREATMENT BALANCE TESTS -- POPULATION MARGIN */
/*---------------------------------------------------------------------------*/

use pop_bw_data, clear

local vars turnout1 effround1 vote_share_t1 area_agr r_total e_total ///
	local_cop margin_local_cop
local n_vars : word count `vars'

matrix pretreat_pop = J(`n_vars',7,.)

forval i = 1/`n_vars' {
	qui rdrobust `: word `i' of `vars'' pop_margin
	storeinv `i' pretreat_pop
	}

matrix rownames pretreat_pop = turnout eff_num_par vote_margin_t1 agr_area ///
	revenues expenditures copartisan copartisan_margin

* show results
matrix list pretreat_pop

	
/*---------------------------------------------------------------------------*/
/* TABLE A7: TESTS OF THE MANIPULATION OF THE RUNNING VARIABLE */
/*---------------------------------------------------------------------------*/

matrix res_mcc = J(6,7,.)

global salary_bw = .186

use electoral_data, clear

* row 1: vote margin
DCdensity vote_margin, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) nograph
	drop Yj-se_
storemcc vote_margin 1 res_mcc

* row 2: vote margin within the balanced population window
DCdensity vote_margin if pop_margin < $salary_bw & pop_margin > -$salary_bw, ///
	breakpoint(0) generate(Xj Yj r0 fhat se_fhat) nograph
	drop Yj-se_
storemcc vote_margin 2 res_mcc

* row 3: vote margin below the population salary threshold
DCdensity vote_margin if pop_margin < 0, breakpoint(0) ///
	generate(Xj Yj r0 fhat se_fhat) nograph
	drop Yj-se_
storemcc vote_margin 3 res_mcc

* row 4: vote margin above the population salary threshold
DCdensity vote_margin if pop_margin > 0, breakpoint(0) ///
	generate(Xj Yj r0 fhat se_fhat) nograph
	drop Yj-se_
storemcc vote_margin 4 res_mcc

* row 5: population margin
DCdensity pop_margin, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) nograph
	drop Yj-se_
storemcc pop_margin 5 res_mcc

* row 6: population margin wihin the balanced population window
DCdensity pop_margin if pop_margin < $salary_bw & pop_margin > -$salary_bw, ///
	breakpoint(0) generate(Xj Yj r0 fhat se_fhat) nograph
	drop Yj-se_
storemcc pop_margin 6 res_mcc

matrix rownames res_mcc = vote_all vote_balanced vote_below vote_above ///
	pop_all pop_balenced
	
* show results
matrix list res_mcc


/*---------------------------------------------------------------------------*/
/* FIGURE A4: BALANCED POPULATION WINDOW */
/*---------------------------------------------------------------------------*/

use pop_bw_data, clear

gen treat = (pop_margin > 0) if pop_margin ~= .

foreach x in sample min_pv w0 {
	gen `x' = .
	}

local w0 500 499 498 497 496 495 494 493 492 491 490 489 488 487 486 485 484 483 482 481 480 479 478 477 476 475 474 473 472 471 470 469 468 467 466 465 464 463 462 461 460 459 458 457 456 455 454 453 452 451 450 449 448 447 446 445 444 443 442 441 440 439 438 437 436 435 434 433 432 431 430 429 428 427 426 425 424 423 422 421 420 419 418 417 416 415 414 413 412 411 410 409 408 407 406 405 404 403 402 401 400 399 398 397 396 395 394 393 392 391 390 389 388 387 386 385 384 383 382 381 380 379 378 377 376 375 374 373 372 371 370 369 368 367 366 365 364 363 362 361 360 359 358 357 356 355 354 353 352 351 350 349 348 347 346 345 344 343 342 341 340 339 338 337 336 335 334 333 332 331 330 329 328 327 326 325 324 323 322 321 320 319 318 317 316 315 314 313 312 311 310 309 308 307 306 305 304 303 302 301 300 299 298 297 296 295 294 293 292 291 290 289 288 287 286 285 284 283 282 281 280 279 278 277 276 275 274 273 272 271 270 269 268 267 266 265 264 263 262 261 260 259 258 257 256 255 254 253 252 251 250 249 248 247 246 245 244 243 242 241 240 239 238 237 236 235 234 233 232 231 230 229 228 227 226 225 224 223 222 221 220 219 218 217 216 215 214 213 212 211 210 209 208 207 206 205 204 203 202 201 200 199 198 197 196 195 194 193 192 191 190 189 188 187 186 185 184 183 182 181 180 179 178 177 176 175 174 173 172 171 170 169 168 167 166 165 164 163 162 161 160 159 158 157 156 155 154 153 152 151 150 149 148 147 146 145 144 143 142 141 140 139 138 137 136 135 134 133 132 131 130 129 128 127 126 125 124 123 122 121 120 119 118 117 116 115 114 113 112 111 110 109 108 107 106 105 104 103 102 101 100 99 98 97 96 95 94 93 92 91 90 89 88 87 86 85 84 83 82 81 80 79 78 77 76 75 74 73 72 71 70 69 68 67 66 65 64 63 62 61 60 59 58 57 56 55 54 53 52 51 50
local nw : word count `w0'
di `nw'

forval i = 1/`nw' {
	local j = `: word `i' of `w0''/1000
	replace w0 = `j' in `i'
	replace sample = .
	replace sample = 1 if pop_margin >= -`j' & pop_margin <= `j'
	foreach x in turnout1 effround1 vote_share_t1 area_agr r_total e_total local_cop ///
		margin_local_cop {
			qui reg `x' treat if sample == 1, vce(cluster year)
			matrix p = r(table)
			local pval = p[4,1]
			gen pv_`x' = `pval' in `i'
			}
	egen min_pv1 = rowmin(pv_*)
	replace min_pv = min_pv1 in `i'
	drop pv_* min_pv1 
	}

* the graph
twoway (scatter min_pv w0, msymbol(Oh) mcolor(black)), scheme(s1mono) ///
	xtitle("Bandwidth Around Salary Threshold (Absolute value)")  ///
	ytitle("Minimum p-value across All Covariates") ///
	aspect(1) yline(.15, lpattern(dash)) xline(.186, lcolor(red))


/*---------------------------------------------------------------------------*/
/* TABLE A8: PROPENSITY TO RUN IN THE NEXT ELECTION */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

matrix running = J(6,7,.)

* row 1: party level, all
qui rdrobust not_run_party vote_margin
storeinv 1 running

* row 2: party level, below salary threshold
qui rdrobust not_run_party vote_margin if pop_margin < 0
storeinv 2 running

* row 3: party level, above salary threshold
qui rdrobust not_run_party vote_margin if pop_margin > 0
storeinv 3 running

* row 4: candidate level, all
qui rdrobust not_run_indiv vote_margin
storeinv 4 running

* row 5: candidate level, below salary threshold
qui rdrobust not_run_indiv vote_margin if pop_margin < 0
storeinv 5 running

* row 6: candidate level, above salary threshold
qui rdrobust not_run_indiv vote_margin if pop_margin > 0
storeinv 6 running

matrix rownames running = party_all party_below_salary party_above_salary ///
	cand_all cand_below_salary cand_above_salary
	
* show results
matrix list running


/*---------------------------------------------------------------------------*/
/* TABLES A9: CENSUS COUNTS MANIPULATION TESTS */
/*---------------------------------------------------------------------------*/

use census_pop_data, clear

reshape long pop_, i(city_code_unq) j(year)
gen pop_margin = (pop_ - 7000)/7000 if pop_ > 3000 & pop_ < 15000

matrix mccs_census = J(3,7,.)

* column 1: 2002 and 2011 census (see Supp. Appendix for explanation)
DCdensity pop_margin, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) ///
	nograph
	drop Yj-se_
storemcc pop_margin 1 mccs_census

* column 2: 2002 census
DCdensity pop_margin if year == 2003, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) ///
	nograph
	drop Yj-se_
storemcc pop_margin 2 mccs_census

* column 3: 2011 census
DCdensity pop_margin if year == 2013, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) ///
	nograph
	drop Yj-se_
storemcc pop_margin 3 mccs_census

matrix rownames mccs_census = 2002_2011_census 2002_census 2011_census

matrix mccs = mccs_census'
* show results
matrix list mccs


/*---------------------------------------------------------------------------*/
/* TABLE A10: WITHIN-LOCALITY ANALYSIS */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

* winner dummy 
gen win_dummy = vote_margin > 0 if vote_margin ~= .
* above salary threshold dummy
gen above = pop_margin > 0 if pop_margin ~= .

gen vint1 = win_dummy*vote_margin
gen pint1 = above*pop_margin

xtset city_code_unq year

xtreg win i.win_dummy##i.above vote_margin vint1 pint1 if vote_margin > -.3 ///
	& vote_margin < .3 & pop_margin > -.186 & pop_margin < .186, r fe

matrix table = r(table)
matrix within = J(4,3,.)
* column 1: below salary threshold
matrix within[1,1] = _b[1.win_dummy]
matrix within[2,1] = _se[1.win_dummy]
matrix within[3,1] = table[4,2]
* column 2: above salary threshold
lincom 1.win_dummy+1.win_dummy#1.above
matrix within[1,2] = r(estimate)
matrix within[2,2] = r(se)
matrix within[3,2] = (1-normal(abs(r(estimate)/r(se))))*2
* column 3: difference
matrix within[1,3] = _b[1.win_dummy#1.above]
matrix within[2,3] = _se[1.win_dummy#1.above]
matrix within[3,3] = table[4,8]
matrix within[4,1] = e(N)

matrix rownames within = Estimate SE p-value N
matrix colnames within = below-salary above-salary difference

* show results
matrix list within


/*---------------------------------------------------------------------------*/
/* TABLE A11: MAYOR'S WEALTH IN 2008 RELATIVE TO LOSERS */
/*---------------------------------------------------------------------------*/

use wealth_data, clear

matrix res_wealth_initv = J(7,3,.)

* column 1: immovable assets
qui rdrobust immovable_2008 vote_margin
stores 1 res_wealth_initv

* column 2: movable assets
qui rdrobust movable_2008 vote_margin
stores 2 res_wealth_initv

* column 3: financial assets
qui rdrobust financial_2008 vote_margin
stores 3 res_wealth_initv

* show results
matrix list res_wealth_initv


/*---------------------------------------------------------------------------*/
/* TABLE A12: MAYOR'S WEALTH IN 2008 BELOW AND ABOVE SALARY THRESHOLD */
/*---------------------------------------------------------------------------*/

matrix res_wealth_initp = J(7,3,.)

* column 1: immovable assets
qui rdrobust immovable_2008 pop_margin if winner == 1
stores 1 res_wealth_initp

* column 2: movable assets
qui rdrobust movable_2008 pop_margin if winner == 1
stores 2 res_wealth_initp

* column 3: financial assets
qui rdrobust financial_2008 pop_margin if winner == 1
stores 3 res_wealth_initp

* show results
matrix list res_wealth_initp


/*---------------------------------------------------------------------------*/
/* TABLE A13: # OF PROCUREMENT CONTRACTS ACROSS SALARY THRESHOLD */
/*---------------------------------------------------------------------------*/

use corruption_data, clear

matrix res_procur_count = J(7,3,.)

* column 1: opaque procedure
qui rdrobust count_opaque pop_margin
stores 1 res_procur_count

* column 2: price per quantity
qui rdrobust count_ppq pop_margin
stores 2 res_procur_count

* column 3: single-bidder contracts
qui rdrobust count_single_bid pop_margin
stores 3 res_procur_count

* show results
matrix list res_procur_count


/*---------------------------------------------------------------------------*/
/* TABLE A14: PROSECUTIONS ACROSS SALARY THRESHOLD */
/*---------------------------------------------------------------------------*/

use prosecution_data, clear

matrix res_pna_ani = J(7,3,.)

* column 1: all prosecutions
qui rdrobust charges pop_margin
stores 1 res_pna_ani

* column 2: ANI charges
qui rdrobust ani_charges pop_margin
stores 2 res_pna_ani

* column 3: # of agencies charged
qui rdrobust org_count pop_margin
stores 3 res_pna_ani

* show results
matrix list res_pna_ani


/*---------------------------------------------------------------------------*/
/* FIGURE A5: SENSITIVITY OF THE INCUMBENCY DISADVANTAGE ESTIMATES 
	ACROSS THE SALARY THRESHOLD */
/*---------------------------------------------------------------------------*/

use electoral_data, clear

matrix multw_below = J(101,3,.)
matrix multw_above = J(101,3,.)

use electoral_data, clear
* for bandwidths [.136, .236] -- i.e. +/- .05 around the balanced pop bandwidth
local steps 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
local nw : word count `steps'
di `nw'

* incumbency disadvantage estimates below and above the salary threshold 
gen bandwidth = .
forval i = 1/`nw' {
	local j = `: word `i' of `steps''/1000
	di `j'
	di `i'
	replace bandwidth = `: word `i' of `steps''/1000 in `i'
	qui rdrobust win vote_margin if pop_margin < 0 & pop_margin > -`j'
	stores_mw `i' multw_below
	qui rdrobust win vote_margin if pop_margin > 0 & pop_margin < `j'
	stores_mw `i' multw_above
	}

* the bootstrapped difference for each bandwidth
gen cond = .
set seed 7
forval k = 1/`nw' {
	local j = `: word `k' of `steps''/1000
	replace cond = .
	replace cond = 0 if pop_margin < 0 & pop_margin > -`j'
	replace cond = 1 if pop_margin > 0 & pop_margin < `j' 
	qui rdrobust win vote_margin if cond == 0
	global h0 = e(h_bw)
	global b0 = e(b_bw)
	qui rdrobust win vote_margin if cond == 1
	global h1 = e(h_bw)
	global b1 = e(b_bw)

	matrix est_diff = J(100,1,.)
	forval i = 1/100 {
		quietly {
			preserve
			bsample, strata(cond)
			rdrobust win vote_margin if cond == 0, h($h0) b($b0)
			local est0 = e(tau_cl)
			rdrobust win vote_margin if cond == 1, h($h1) b($b1)
			local est1 = e(tau_cl)
			matrix est_diff[`i',1] = `est1' - `est0'
			restore
			}
		}
	svmat est_diff
	rename est_diff1 est_diff_`k'
	di `k'
	}

gen diff_l = .
gen diff = .
gen diff_u = .
forval i = 1/101 {
	_pctile est_diff_`i', p(2.5 50 97.5)
	replace diff_l = r(r1) in `i'
	replace diff = r(r2) in `i'
	replace diff_u = r(r3) in `i'
	}
	
*graph the difference
twoway (line diff bandwidth, lcolor(black)) ///
	(line diff_l bandwidth, lpattern(dash) lcolor(black)) ///
	(line diff_u bandwidth, lpattern(dash) lcolor(black)), ///
	scheme(s1mono) ylabel(, nogrid) yline(0, lcolor(gs13)) legend(off) ///
	aspect(0) xtitle(Bandwidth Size Around Salary Threshold) ///
	ytitle("Difference in Incumbency Advantage" "Across Salary Threshold") ///
	yscale(range(-1(.2)1)) ylabel(#7) name(diff_graph, replace) ///
	xline($salary_bw, lcolor(red))
	
* and show the individual estimates
svmat multw_below
svmat multw_above
* calculate the 95% CI
foreach x in below above {
	gen multw_`x'_l = multw_`x'1 - 1.96*multw_`x'2
	gen multw_`x'_u = multw_`x'1 + 1.96*multw_`x'2
	}

* graph estimates below the salary threshold
sum multw_below3
local nmin = r(min)
local nmax = r(max)
twoway (line multw_below1 bandwidth, lcolor(black)) ///
	(line multw_below_l bandwidth, lpattern(dash) lcolor(black)) ///
	(line multw_below_u bandwidth, lpattern(dash) lcolor(black)), ///
	scheme(s1mono) ylabel(, nogrid) yline(0, lcolor(gs13)) legend(off) ///
	aspect(0) xtitle(Bandwidth Size Below Salary Threshold) ///
	ytitle("Incumbency Advantage Below" "Salary Threshold") ///
	yscale(range(-1(.2)1)) ylabel(#7) name(below_graph, replace) ///
	xline($salary_bw, lcolor(red)) text(.95 .225 "Min N: `nmin'" .8 .225 "Max N: `nmax'")

* graph estimates above the salary threshold
sum multw_above3
local nmin = r(min)
local nmax = r(max)
twoway (line multw_above1 bandwidth, lcolor(black)) ///
	(line multw_above_l bandwidth, lpattern(dash) lcolor(black)) ///
	(line multw_above_u bandwidth, lpattern(dash) lcolor(black)), ///
	scheme(s1mono) ylabel(, nogrid) yline(0, lcolor(gs13)) legend(off) ///
	aspect(0) xtitle(Bandwidth Size Above Salary Threshold) ///
	ytitle("Incumbency Advantage Above" "Salary Threshold") ///
	yscale(range(-1(.2)1)) ylabel(#7) name(above_graph, replace) ///
	xline($salary_bw, lcolor(red)) text(-.8 .225 "Min N: `nmin'" -.95 .225 "Max N: `nmax'")
	
gr combine below_graph above_graph diff_graph, scheme(s1mono)

graph drop _all
